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ABSTRACT 

We present the first demonstration of the negative effective magnetic pressure instability in direct numerical 
simulations of stably stratified, externally forced, isothermal hydromagnetic turbulence in the regime of large 
plasma beta. By the action of this instability, an initially uniform horizontal magnetic field forms flux concen- 
trations whose scale is large compared to the turbulent scale. We further show that the magnetic energy of these 
large-scale structures is only weakly dependent on the magnetic Reynolds number, provided its value is large 
enough for the instability to be excited. Our results support earlier mean-field calculations and analytic work 
which identified this instability. Applications to the formation of active regions in the Sun are discussed. 
Subject headings: magnetohydrodynamics (MHD) - Sun: dynamo - sunspots - turbulence 
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1. INTRODUCTION 

The solar convection zone is highly turbulent and mixing is 
expected to be efficient. Nevertheless, the Sun displays co- 
herent structures encompassing many turbulent eddy scales. 
A well-known example is the large-scale magnetic field of 
the Sun that is antis ymmetric about the eq uator and shows 
a 22 year solar cycle ( Stenflo & Vogell [l986). Another promi- 
nent example in the Sun is the emergence of active regions. 
It is generally believed that active regions are the result of 
some non-axisymmetri c instability of ~ 100 kG magnetic 
fields in the tachocline (iGilman & Dikpati 2000; iCallv et alj 
120031: iParfrev & Menoul 120071) . Ho wever, the existenc e of 
such strong fields remains debatable (Brandenburg 20051). 

A powerful tool for understanding the emergence of 
such large-scale structures fr om a turbulent background 
is mean-field dynamo theory dMoffattl 119781 ; iParke 
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iKrause & Radlerll980l) . With the advent of powerful comput- 
ers and numerical simulation tools, it has become possible to 
confront many of the mean-field predictions with direct sim- 
ulations (Brandenburg & Subramanian 2005). Here we con- 
sider the idea that statistically steady, stratified, hydromag- 
netic turbulence with an initially uniform magnetic field is 
unstable to the negative effective magnetic pressure instabil- 
ity (NEMPI). This instability is caused by the suppression of 
turbulent hydromagnetic pressure (the isotropic part of com- 
bined Reynolds and Maxw e ll stresses) by the mean magneti c 
field dKleeorinet alj fl990l: iRogachevskii & Kleeoriij|2007T) . 
At large Reynolds numbers and for sub-equipartition mag- 
netic fields, the negative turbulent contribution can become 
so large that the effective mean magnetic pressure (the sum of 
turbulent and non-turbulent contributions) appears negative. 
In a stratified medium, this results in the excitation of NEMPI 
that causes formation of large-scale inhomogeneous magnetic 
structures. NEMPI is similar to the large-scale dynamo insta- 
bility, except that it only redistributes the total magnetic flux, 
creating large-scale concentrated magnetic flux regions at the 
expense of turbulent kinetic energy. 

Historically, the mag netic suppress ion of the Reynold s 
stress was first found bvlRadlerl ( 119741) and iRiidigerl (119741) . 
Later. iRtidiger et aT] d 1 986b considered the Maxwell stress and 
found the mean effective magnetic tension to be suppressed 



by mean fields. However, these calculations were based on 
quasi-linear theory which is only valid at low fluid a nd mag- 
netic Reynolds numbers. Kleeorin et al.l (119901119961) consid- 
ered the combined Reynolds and Maxwell stresses at large 
Reynolds numbers and found a sign reversal of the effective 
mean magnetic pressure. This result is based on the r approx- 
imation, an d has been corroborated using th e renormalization 
procedure (Kle eorin & Ro gachevskii 1994]). 

The magnetic suppression of the combined Reynolds and 
Maxwell stresses is quantified in terms of new turbulent 
mean-field coefficients that relate the components of the sum 
of Reynolds and Maxwell stresses to the mean magnetic 
field. These coefficients depend on the magnetic field and 
have now been determined in direct numerical simulations 
(DNS) for a broad range of dif ferent cases, in cluding unstrati- 
fied forced turbulence ( Brandenbur g et alJ2010l). isothe rm ally 
stratified forced turbulence ( Brandenburg e t all 1201 11 here - 
after BKKR), and turbulent convection (Kapvla et alj|20llb . 
These simulations have clearly demonstrated that the mean 
effective magnetic pressure is negative for magnetic field 
strengths below about half the equipartition field strength. 
However, these DNS studies did not find the actual instability. 

With a quantitative parameterization in place, it became 
possible to build mean-field models of stratified turbulence 
which clearly demonstrate exponential growth and satura- 
tion of NEMPI. In view of applications to the formation of 
active regions in the Sun, such simulatio ns were originally 
done for an adiabatically stratified layer (Brandenbu rg et all 
2010). In addition, mean-field studies showed the exis- 
tence of NEMPI e ven for isothermal stably stratified layers 
dKemel et al.l201 lb . This last result turned out to be important 
because it paved the way for this Letter where we demonstrate 
NEMPI through DNS. Once we establish the physical reality 
of this effect, it would be important to apply it to realistic solar 
models which include proper boundary conditions, realistic 
stratification, convective flux, and radiation transport. How- 
ever, at this stage it is essential to isolate NEMPI as a physical 
effect under conditions that are as simple as possible. 



2. THE MODEL 

We consider a domain of size L^x L 



x L z in Cartesian co- 
ordinates, (x, y, z), with periodic boundary conditions in the 
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AB y / B C qo in the xz plane for Rejvf = 6 and Bo/B cq o = 0.05, showing a descending "potato sack" structure. Time is in units of rto (lower 



x and y directions and stress-free perfectly conducting bound- 
aries at top and bottom (z — ±L z /2). The volume-averaged 
density is therefore constant in time and equal to its initial 
value, po = (p) . We solve the equations for the velocity U, 
the magnetic vector potential A, and the density p, 



p^- = J x B c 2 Vp- 
H Dt 3 ' 



V-(2upS)+p(f+g), (1) 



| = t/xB + „V 2 A, 



dp 
dt 



-V-pU, 



(2) 



(3) 



where v is kinematic viscosity, r\ is magnetic diffusivity, B = 
B a + V x A is the magnetic field, B a = (0, B , 0) is the im- 
posed uniform field, J = V x B/po is the current density, po 
is the vacuum permeability, S,j = ^(Uij +Uj t i)— | <5jjV-Z7 
is the traceless rate of strain tensor, and commas denote par- 
tial differentiation. The forcing function / consists of ran- 
dom, white-in-time, plane non-polarized waves with an aver- 
age wavenumber kf = 15 fci, where k\ — 2ir/L z is the lowest 
wavenumber in the domain. The forcing strength is such that 
the turbulent rms velocity is approximately independent of z 
with u rms = (u, 2 ) 1 / 2 R3 0.1 c s . The gravitational acceleration 
g = (0, 0, —g) is chosen such that k\H p = 1, which leads to 
a density contrast between bottom and top of exp(27r) s» 535. 
Here, H p = c^/g is the density scale height. 

Our simulations are characterized by the fluid Reynolds 
number Re = u Tms /vkf, the magnetic Prandtl number PrM = 
v/rj and the magnetic R eynolds number ReM = RePi"M- 
Following earlier work (Brandenb urg et ai]|201 ll) . we choose 
Pi"m = 0.5 and Re^f in the range 0.7—74. The magnetic field 
is expressed in units of the local equipartition field strength 
near the top, B cq = y/popu lnls , while Bq is specified in 
units of the averaged value, £> cq0 = ^poPo u lms . We moni- 
tor AB y = By — Bq, where B y is an average over y and a 
certain time interval At. Time is expressed in eddy turnover 
times, Tto = (urmsfc/) -1 ' Occasionally, we also consider 
the turbulent-diffusive timescale, rtd = (vtokf)^ 1 , where 
? /to = w rms/3fc/ is the estimated turbulent magnetic diffu- 
sivity. Another diagnostic quantity is the rms magnetic field 
in the k = ki Fourier mode, B%, which is here taken as an 
average over 2 < kiz < 3, and is close to the top at k\Z = n. 
(Note that B\ does not include the imposed field Bq at k — 0.) 

The simulations are performed with the PENCIL Code, 1 

1 |http : / /pencil-code, googlecode . com| 
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FIG. 2. — Dependence of Bi/B cq (averaged over 2 < k\z < 3)onReM 
for Bo/BoqO = 0.05. 

which uses sixth-order explicit finite differences in space and 
a third-order accurate time stepping method. We use numer- 
ical resolutions of 128 3 and 256 3 mesh points when L x — 
L y = L z , and 1024 x 128 2 when L x = 8L y = 8L 2 . To 
capture mean-field effects on the slower turbulent-diffusive 
timescale, which is T t d/T to = 3fc 2 //c 2 times slower than 
the dynamical timescale, we perform simulations for several 
thousand turnover times. 

3. RESULTS 

The NEMPI phenomenon is already quite pronounced at 
intermediate values of Re^/ <; 3; see Figure [T] where we 
show By averaged over y and At w 80r to (denoted by B y ) at 
selected times during the first 2000 turnover times. The ReM 
dependence of Bi/B cq is shown in Figure [2] For Re m > 3, 
the Re^f dependence is relatively weak, which suggests that 
NEMPI is a genuine high Reynolds number effect. 

The results for the case of Re^j = 6 show strong similar- 
ities to earlier mean-field simulations. During the first 500 
turnover times, flux concentrations form first near the surface, 
but at later times the location of the peak magnetic field moves 
gradually downward. This phenomenon is a direct conse- 
quence of the negative effective magnetic pressure, making 
such structures heavier than their surroundings. Their shape 
resembles that of a falling "potato sack" and has been seen in 
numerous mean-field calculations during the nonlinear stage 
of NEMPI dBrandenbure et akllloTotlkapvla et al.ll20TH) . Us- 
ing the technique described in BKKR, we have found that for 
Re jt/ 1.1 and Pim = 0.5, the effective magnetic pressure 
has a negative minimum. 
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FIG. 3.— AB H /BeqO in the xz plane for Re M = 36 and B /B eq0 = 0.02, 0.05, 0.09, and 0.18 
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FIG. 4. — Bifurcation diagram showing Bi/B cq vs. Bq/B cci q (and vs. 
So / B cq on the upper abscissa) for ReM = 36. 

As expected from theory and mean-field calculations, 
NEMPI is only excited in a certain range of field strengths. 
In particular, only for B /B cq( ) between 0.02 and 0.2 do 
we see large-scale magnetic structures. This is shown in 
Figure [3] where we see B y , again averaged over y and a 
time interval At « 800r to , in which the field is statistically 
steady. The clearest flux structure formations are seen for 
Bo/B cq o m 0.05. However, even for this case the flux con- 
centrations are barely visible in a single snapshot. This has 
been one of the reasons why NEMPI has not been noticed 
before in DNS. An additional handicap was that the simula- 
tions of BKKR used a smaller scale separation ratio of only 5, 
which is nevertheless still sufficient for determining the gov- 
erning mean-field coefficients and allows one to reach larger 
values of ReM- 

In Figure 2] we plot Bi/B cq as a function of Bo/B cq o, 
showing a peak at Bo/B cq o w 0.05. We recall that B cq ap- 
plies here to the location 2 < k\z < 3 where B\ has been 
evaluated, and there we have B eq / B eq o ~ 0.3. The fact 
that large-scale flux concentrations develop only for a cer- 
tain range of imposed field strengths supports our interpre- 
tation that they are caused by NEMPI and not, for example, 
by some yet unknown dynamo mechanism. In all these cases, 
B — Bq grows rapidly and reaches a saturation field strength 
that is independent of Bq, provided Rejv/ > 35. This sug- 
gests that this field is produced by small-scale dynamo ac- 
tion and not just by field line tangling. Another piece of evi- 
dence of the physical reality of NEMPI is shown in Figure |5J 
where we see that Bi does indeed increase exponentially for 
the first 2000 turnover times, corresponding to about 3rtd- 
The growth rate is w OArjtokf, which is much less than t,." 1 , 



Fig. 5. — Slow exponential growth of the mean magnetic field in 2 < 
k\z < 3 for Bo/B cq o = 0.05, corresponding to Bo/B cq = 0.16, for 
Rejvf = 36 (solid line). The growth rate is f» OA'rjtok'i (dotted line). Note 
that the total rms field, B Tms / B cq (dotted line) saturates much faster, as 
indicated by to upper abscissa. 

but entirely c ompatible with mean-field calculations (BKKR; 
Kapyla et al.ll201lb . 

Finally, to investigate the effects of the domain aspect ratio 
on the instability, we perform a calculation with Bq /B eq o = 
0.05, ReM = 36, and change L x to 167r/fci. We find that the 
most unstable mode has a wavelength approximately equal 
to L z s» 6H p . This result i s also in agreement with mean- 
field models (e.g., Fig. 14 of iKapyla et alJl201lh . The large- 
scale flux concentrations have an amplitude of only w 0.1£? cq 
and are therefore not seen in single snapshots, where the field 
reaches peak strengths comparable to B cq . Furthermore, as 
for any linear instability, the flux concentrations form a repet- 
itive pattern, but this might be an artifact of idealized condi- 
tions. 

4. CONCLUSIONS AND DISCUSSION 

The present simulations have, for the first time, demon- 
strated conclusively that NEMPI can operate in hydromag- 
netic turbulence under proper conditions, namely, strong strat- 
ification, sufficient scale separation (here fc//fci = 15), and a 
mean field in an optimal range (here w 0.15£? oq ; see Fig- 
ure|l]i. This instability has so far only been seen in mean-field 
simulations. By contrast, the present simulations are com- 
pletely free of any mean-field consideration. 

The instability is argued to be a consequence of the reduc- 
tion of turbulent hydromagnetic pressur e by a mean magneti c 
field and can be understood as follows (Kle eorin et alj 1996). 
The combined Reynolds and Maxwell stress is pUiUj — 
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5ijb 2 /2fi , where u and b are velocity and mag- 
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netic fluctuations, respectively, and overbars indicate averag- 
ing. For isotropic turbulence, the turbulent hydromagnetic 
pressure is then Pturb = \{pu 2 + b 2 /2/xq). On the other 

hand, the total turbulent energy Etmb = \{p u2 + b 2 /fio), 
is nearly conserved because a uniform mean magnetic fi eld 
does not perform any work; see Brandenbur g et alJ (120101) for 
a numerical demonstration. The presence of an additional 
1/2 factor in front of the b 2 //io term in the expression for 
-Pturb, but not in that for Pturb, implies that the generation 
of magnetic fluctuations results in a reduction of Pturb, i-e- 
Pturb = |(2Pturb ^ b 2 /2/io)- For anisotropic turbulence this 
negative contribution becomes larger. This physical effect is 
independent of stratification, but to obtain an instability one 
needs strong stratification. 

We speculate that, in the solar context, NEMPI plays a role 
in formation of active regions from mean fields generated by 
the solar dynamo. Let us now ask whether this instability 
alone can describe the formation of active regions at the so- 
lar surface. Clearly, the flux concentrations we observe are 
not strong enough to be noticeable without averaging, while 
the active regions in the Sun are seen without averaging. This 
suggests that there may be additional mechanisms at work. 
One possibility is that of the magnetic suppression of the con- 
vective heat flux that has been invoked to ex plain the forma- 
tion of sunspots (Kitchatinov & Mazur 2000). 

When the mean magnetic field becomes larger than the 
equipartition field strength, i.e., when NEMPI does not work, 
and the characteristic spatial scale of the magnetic field is 
smaller than the density h eight, it is ins tead the Parker mag- 
netic buoyancy instability (Parker 1966) that is excited. 

The presence of a vertical field might also have a strong 
effect. Indeed simulations of convection with an imposed ver- 
tical field have produ ced a segregation into magnetized and 
unmagnetized regions dTaoetal.ll 99 8; Kitiashvili et al. 2010) 
with formation of flux concentrations strong enough to be no- 
ticeable even without averaging. On the other hand NEMPI 
might become more powerful at stronger stratification. In- 
creased strat ification clearly ha s an enhancing effect on the 
growth rate (Kernel et al.H2.Ql ll) . but the effect on the satura- 
tion level has not yet been quantified. Furthermore, the inter- 
play between NEMPI and the other effects also needs to be 
investigated. 

Our work has established a close link between what can 
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be expected from mean-field studies and what actually hap- 
pens in DNS. This correspondence is particularly important 
because DNS cannot reach solar parameters in any conceiv- 
able future. Hence a deeper understanding of solar convection 
can only emerge by studying mean-field models on the one 
hand and to determine turbulent mean-field coefficients from 
simulations on the other hand. This concerns not only the de- 
pendence of the mean-field coefficients on parameters such 
as magnetic Reynolds and Prandtl numbers and scale sepa- 
ration ratio, but also the details of the source of turbulence. 
In particular, it has already been shown that the negative ef- 
fective magnetic pressure effect is not unique to f orced turbu- 
lence, but it also occurs in turbulent convection Kapvla et al. 
(1201 11) . and thus in an unstably stratified layer. We empha- 
size that the present work demonstrates the predictive power 
of mean-fiel d theory at an advan ced level where quasi-linear 
theory fails (Rudiger et al.jl2.01 ll) but the spectral r approach 
(Rogachevskii & Kleeorin 2007) has proven useful. 

More work using mean-field models is needed to eluci- 
date details of the mechanism of NEMPI. For example, naive 
thinking suggests that the onset of NEMPI should occur at 
the depth where the effective magnetic pressure is minimum, 
but both mean-field models and DNS show that this is not 
the case. At least at early times, NEMPI appears most pro- 
nounced at the top of the domain, while the effective mag- 
netic pressure is usually most negative at the bottom. On the 
other hand, the instability is a global one and local consid- 
erations such as these are not always meaningful. Another 
question is what happens when the imposed field is replaced 
by a dynamo-generated one. In that case, the turbulence may 
be helical and new terms involving current density can occur 
in the expression for the mean-field stress. Again, such possi- 
bilities are best studied using first the mean-field approach. 
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